\defmodule {NI3}
This class implements the algorithm NI3 (protected Newton-Raphson method).
  The root is found with accuracy tolerance.

\bigskip\hrule

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        NI3
 * Description:
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Université de Montréal
 * Organization: DIRO, Université de Montréal
 * @author       Nabil Channouf
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.probdistmulti.norta;
\begin{hide}
import umontreal.iro.lecuyer.probdist.*;
\end{hide}

public class NI3 extends NortaInitDisc \begin{hide}
{
   private double tolerance; /* Desired accuracy for the root-finder
   			        algorithm (epsilon in paragraph
   			        "Method NI3" of section 3 in paper).*/
\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}

\begin{code}

   public NI3 (double rX, DiscreteDistributionInt dist1,
               DiscreteDistributionInt dist2, double tr, double tolerance)\begin{hide}
   {
      super(rX, dist1, dist2, tr);
      this.tolerance = tolerance;
      computeParams();
   }\end{hide}
\end{code}
\begin{tabb}
   Constructor of the class NI3 with the target rank correlation rX,
       the two discrete marginals dist1 and dist2,
     the parameter for truncation tr and the specific parameter
      tolerance which corresponds to epsilon in the paper
      (paragraph "Method NI3" of section 3).
\end{tabb}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}
\begin{code}

   public double computeCorr() \begin{hide}
   {
      final double ITMAX = 100; // Maximum number of iterations.
      double xl, xh;  /* Left and right endpoints of the root bracket at
			 all iterations. */
      double b = 0.0; // The returned solution.
      double f, df;   // Function and its derivative evaluations.
      double dx;      // Correction term.
      /** f, df, dx correspond to f(rho_k), f'(rho_k) and f(rho_k)/f'(rho_k),
          respectively, in the paper (paragraph "Method NI3"of section 3). */
      double dxold; // The root correction at one iteration before the last one.
      double temp;// The root at one iteration before the last one.
      double ccc = rX * sd1 * sd2 + mu1 * mu2; // Precompute constant.

      if (rX == 0.0)
         return 0.0;
      if (rX > 0.0) {              // Orient the search
         xl = 0.0;
         xh = 1.0;
      } else {
         xl = -1.0;
         xh = 0.0;
      }

      b = 2 * Math.sin (Math.PI * rX / 6); // Initial guess
      dxold = xh - xl;
      dx = dxold;
      f = integ (b) - ccc;
      df = deriv (b);

      for (int i = 1; i <= ITMAX; i++) { // Begin the search
         if ((((b - xh) * df - f) * ((b - xl) * df - f) > 0.0)
               || (Math.abs (2.0 * f) > Math.abs (dxold * df))) {
            // Do bisection if solution is out of range
            // or not decreasing fast enough
            dxold = dx;
            dx = 0.5 * (xh - xl);
            b = xl + dx;
            if (xl == b)      // Change in root is negligible.
               return b;      // Accept this root
         } else {
            dxold = dx;
            dx = f / df;
            temp = b;
            b -= dx;
            if (temp == b)
               return b;
         }
         if (Math.abs (dx) < tolerance)
            return b;          // Convergence check
         f = integ (b) - ccc;
         df = deriv (b);
         if (f < 0.0)
            xl = b;            // Maintain the brackets on the root
         else
            xh = b;
      }
      return b;
   }\end{hide}
\end{code}
\begin{tabb} Computes and returns the correlation $\rho_Z$ using the algorithm NI3.
\end{tabb}
\begin{code}

   public String toString()\begin{hide}
   {
    // To display the inputs.
      String desc = super.toString();
      desc += "tolerance : " + tolerance + "\n";
      return desc;
   }\end{hide}
\end{code}

\begin{code}\begin{hide}
}\end{hide}
\end{code}
